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ABSTRACT 


We study the interaction of feedback from active galactic nuclei (AGN) and a multi-phase interstel¬ 
lar medium (ISM), in simulations including explicit stellar feedback, multi-phase cooling, accretion-disk 
winds, and Compton heating. We examine radii ^ 0.1 — lOOpc around a black hole (BH), where the accre¬ 
tion rate onto the BH is determined and where AGN-powered winds and radiation couple to the ISM. We 
conclude: (1) The BH accretion rate is determined by exchange of angular momentum between gas and 
stars in gravitational instabilities. This produces accretion rates ^ 0.03 — lMQyr“\ sufficient to power 
luminous AGN. (2) The gas disk in the galactic nucleus undergoes an initial burst of star formation fol¬ 
lowed by several Myrs where stellar feedback suppresses the star formation rate (SFR). (3) AGN winds 
injected at small radii with momentum fiuxes ^ Lagn/c couple efficiently to the ISM and have dramatic 
effects on ISM properties within ~ 100 pc. AGN winds suppress the nuclear SFR by factors ^10 — 30 and 
BH accretion rate by factors ^3 — 30. They increase the outfiow rate from the nucleus by factors ^10, 
consistent with observational evidence for galaxy-scale AGN-driven outflows. (4) With AGN feedback, 
the predicted column density distribution to the BH is consistent with observations. Absent AGN feed¬ 
back, the BH is isotropically obscured and there are not enough optically-thin sightlines to explain Type-I 
AGN. A ‘torus-like’ geometry arises self-consistently as AGN feedback evacuates gas in polar regions. 

Key words: galaxies: formation — galaxies: evolution — galaxies: active — star formation: general — 
cosmology: theory 


1 INTRODUCTION 


The masses of super-massive black holes (BHs) correlate with var- 


ious host galaxy bulge properties ([Magorrian et al.|1998|[Ferrarese 

& Merritt|2000| 

[Gebhardt et al.[2000| [Hopkins et al.[2007|[Aller & 

Richstone|2007[ 

, Feoli & Mancini[2009| Kormendy et al.|201 1| for 


relations (relative to other galaxy properties; [Hopkins et al.|2009b[ 
[Kormendy & Ho|r2013| l, together with constraints indicating that 
most BH mass is assembled in an optically bright quasar phase 
(|Soltan|1982||Salucci et al.|1999||Yu & Tremaine |2002||Hopkins| 

let al.|2006bK has led to the development of models where large- 

scale effects of feedback from accretion self-regulate BH growth 
at a c ritical mass ([Silk & Rees|19^ |King|200^ jPi Matteo et al.j 
|2005[|Murray et al.|2005| ). Gas inflows triggered by some process 

fuel rapid BH growth in a nuclear starburst ( jDiamond-Stanic &| 
|Rieke|20T^ [Mushotzky et al.|2014j l, until feedback begins to ex¬ 
pel nearby gas and dust. This “blowout” results in a short-lived, 
bright optical quasar that, having expelled its fuel supply, fades and 
leaves a remnant on the observed BH-host correlations ( [Hopkins j 
[et al.[2005a[c| ). This general scenario has been able to explain many 
quasar observables, including luminosity functions, lifetimes, and 


BH mass functions ([Hopkins et al. 

|[2005b[[2006c[ [2008b 2009a 

[Volonteri et al. [[20061 [Menci et al.[ 

[20031 ISomerville et al.||2008 
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[Lapi et al.|2006[ [Tortora et al.|2009| l. It has also been speculated 
that this feedback might ultimately have a large impact through¬ 
out the AGN host galaxy, expelling or heating gas and explaining 


the rapid quenching of star formation in massive galaxies ([Granato 

et al.| 

2004| [Scannapieco & Oh|2004| [Croton et al.|2006l [Hopkins 

et al.| 

[2008a| Antonuccio-Delogu & Silk[[2008[l, and considerable 


observational evidence has emerged for this in recent years (see 


e.g.[Nesvadba et al. 

2010| Cimatti et al.[2013|[LaMassa et al.[2013| 

[Shimizu et al.[2015 

IGuillard et al.|2015||Alatalo et al.|2015li. 


High-velocity outflows can be driven from the BH accretion 
disk by a variety of physical processes including, e.g., radiation 
pressure on lines and dust, magnetic processes, or Compton heat- 


ing (see e.g. [Blandford & Payne[[1982[ [Begelman[ 

19851 [Chang 

et al.[[1987 

[Sanders et al.[[1988| [Konigl & KartjeJ] 

994| [Murray 

etal.|1995l 

Elvis [2000||Proga 

2000 |2007| Silk|2005| 

Murray et al. 

20051 [Bate] 

leldor et al. 120071 

Tortora et al.|2009[l. These manifest 


themselves observationally as ultra-fast outflows ( [Tombesi et al.[ 
[2010[[20T3[[2015[ e.g.), the broad emission line regions and broad 


absorption line quasars (e.g. Weymann et al. 

1981|[de Kool et al.[ 

|200Tl|Gabel et al.|2006 |Ganguly et al.|2007 

j-n 1 rvi 1 1 —1\ : 

, more moderate ve- 


locity outflows (v ~ 10"^ — 10"^ km s ^) associated with the narrow 


line region ([Laor et al.[[1997[ [Crenshaw et al.[[2000[ [Steenbrugge[ 

[et al.[2005|[Krongold et al.[2007[l, as well as quasar ab sorption and 

occultation systems (e.g. |McKernan & Yaqoob| 1998 [[Turner et al.] 

[2008[[Miller et al.|2008[ >. Observations on galaxy scales have also 

provided strong evidence for powerful molecular, atomic, and ion¬ 
ized outflows with velocities ~ 1 — 5 x 10^kms“\ outflow rates 
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2 Hopkins et al. 


up to ~ 100 — 1000 times the BH accretion rate, and spatial ex- 


nawi|20091jMoe et al.|20091jWild et al.|20091 [Fischer et al.|20101 

Humphrey et al.|2010[jDunn et al.|2010[ [Bautista et al.|2010| 

iFer- 

uglio et al.| 

[2010|[Sturm et al.|2011| 

Rupke & Veilleux|2011|| 

Coil 

et al.|2011 

Greene et al.|2011[[2012 

[Faucher-Giguere et al.|2012| 


[ et al.|2014||2015[|Zakamska & Greene|2014[|Zakamska et al. 


2015| ). In some cases, however, it remains unclear to what extent 
these outflows are driven by AGN activity vs. star formation. 

The physics of how AGN-powered outflows interact with the 
ISM and affect the fueling of the AGN itself - how inflow and out¬ 
flow are governed on scales between the small-scale viscous accre¬ 
tion disk (<^ 0.1 pc) and the galaxy proper (> 0.1 kpc) - remains 
highly uncertain. There have been many theoretical studies of dif¬ 
ferent “modes” of AGN feedback (see references above, as well as 
Ostriker et al.|2010[|Ciotti et al.|2010[[Choi et al.|2014||2015[|Stein-| 

born et aLl|2015[ and references therein). However, to date, most 
of these studies have treated the interstellar medium with relatively 
simple “sub-grid” prescriptions that ignore the additional complica¬ 
tions introduced by stellar feedback, interstellar turbulence, and/or 
the small-scale phase structure of the medium around an AGN. In 
order to build on these models and model the interaction of AGN 
outflows and the ISM with greater fldelity, it is critical to include 
both a realistic description of the physics of the ISM, star forma¬ 
tion, and stellar feedback, as well as a plausible description of AGN 
feedback mechanisms. 

Towards this end, in this paper we use a suite of numerical 
simulations to study the interaction of quasar-driven winds and a 
multi-phase ISM. In a series of papers ( [Hopkins et al.|201 l]|20I2d| ) 
(hereafter Papers I & II, respectively), we have developed a new 
set of numerical methods to explicitly model some of the key pro¬ 
cesses that shape the multi-phase ISM; the simulations include 
physically motivated, but still subgrid, treatments of stellar radi¬ 
ation pressure, HII photoionization and photoelectric heating, and 
the heating, momentum, and mass deposition by supernovae (SNe) 
and stellar winds. The feedback is tied to the young stars with en¬ 
ergetics and time-dependence taken directly from stellar evolution 
models - this is particularly important in galactic nuclei, since the 
dynamical times become shorter than stellar evolution timescales. 


In a series of papers ([Hopkins et al.[[2014| [Muratov et al.j 

[2015 

Onorbe et al.[[2015| [van de Voort et al.[[2015| [Ma et al.j 

2015 

201 6| [Faucher-Giguere et al.[[2015[l, we showed that, on galactic 


scales, these models produce a quasi-steady ISM in which molecu¬ 
lar clouds form and disperse rapidly, with phase structure, turbulent 
velocity dispersions, and disk and GMC properties in reasonable 
agreement with observations. Here, we combine these models with 
models for AGN accretion and feedback via both Compton heating 
and high-velocity winds from the AGN accretion disk, and examine 
how various forms of AGN feedback affect black hole accretion, 
AGN obscuration, and the generation of galaxy-scale outflows. We 
focus on scales of ~ 0.1 — lOOpc, where the accretion rate onto 
the black hole is determined, and where AGN-powered winds and 
Compton heating couple to the ISM. For comparison, the “radius 
of influence” of the BH, inside of which it dominates the potential, 
is ~ GMbh/ct^ ~ 5 pc in the case we study. This is the first in a new 
series of papers so we highlight a few of the key results but leave 
more detailed studies for future work. 

The remainder of this paper is organized as follows. 0 sum¬ 
marizes our galaxy models and our treatment of radiative cooling, 
star formation and BH growth, and stellar and AGN feedback. g 
summarizes the results of simulations with stellar feedback only. 


Table 1. Simulations 


Model 


r\E 

fi 

vbal 

Notes 

no_BAL 

0 

0 

0 

0 

no AGN FB 

v5000 

1 

0.008 

6.0 

5,000 

“default” 

v5000_hiP 

10 

0.08 

60 

5,000 

high-momentum 

vSOOOJoP 

0.1 

0.0008 

0.6 

5,000 

low-momentum 

V30000 

1 

0.05 

1.0 

30,000 

high-energy 

v500 

1 

0.0008 

60 

500 

low-energy 

v5000_C 

1 

0.008 

6.0 

5,000 

-i-Compton heating 

v5000 iso 

1 

0.008 

6.0 

5,000 

isotropic winds 


Parameters describing the simulations in the text: Each employs a 
gas particle mass of 13.5 h~^ Mq and minimum SPH smoothing 
length of 0.0014/?“^ pc. Additional simulations for numerical tests 
are in Appendix [C] 

(1) Model name 

(2) r]p: Momentum-loading of BAL wind feedback {p = r]pL/c) 

(3) pe: Energy-loading of BAL wind feedback (E = tjeL) 

(4) /3: Mass-loading (3 = Mbal/4^bh (determined by pp & pe) 

(5) vbal- AGN wind launching velocity at the simulation 
resolution (in kms“^; determined by pp Slije) 


while ^ compares these results to simulations that include AGN 
feedback. ^ summarizes and discusses our key results. A series 
of Appendices contain key technical results. Appendix A describes 
our implementation of BH feedback. Appendix B summarizes the 
effects of including short timescale variability in the assumed BH 
accretion rate. Appendix C describes convergence tests and the ef¬ 
fects of using alternate numerical methods. Appendix D shows that 
in-shock cooling does not compromise our results. 


2 THE SIMULATIONS 

The simulations were performed using the GIZMO code ( jHopkinsj 
|2015j l. GIZMO is a multi-method code which can be run with any of 
several hydro solvers; here we run the code in its smoothed-particle 
hydrodynamics (SPH) mode, speciflcally in the “pressure-entropy” 
(“P-SPH”) form which includes several improvements relative to 
older SPH implementations. Speciflcally, this is a heavily mod- 
ifled version of the parallel TreeSPH code GADGET-S ( jSpringel 
20051, in a fully conservative formulation ( jSpringel & Hernquist 
20021 which is also density-independent in a manner that allows 
contact discontinuities and improved fluid mixing ( |Hopkins|2013| 
see Appendix]^. The artiflcial viscosity, adaptive timestepping, 
and smoothing kernel are updated following jHopkinsj ( |20I3| ). The 
galaxy models and the treatment of star formation and stellar feed¬ 
back are described in detail in Paper I (Sec. 2 & Tables 1-3) and 
Paper II (Sec. 2). We briefly summarize the salient properties here. 

2.1 Initial Conditions 


The initial conditions are a gas-rich nuclear disk in a massive 
galaxy, drawn f rom the large parameter survey of [Hopkins 


Quataert {201 Oa I .We consider a BH (initial Mbh = 3 x 10' Mq) in 


Hernquist (19901 stellar bulge (Mbuige = IO^^Mq, isotropic orbits 


and scale-length a = 1.7 kpc) and halo (Mhaio = 2 x 10 Mq, with 
virial radius, concentration, and velocity appropriate at z = 0). The 
BH is surrounded by an exponential nuclear disk of gas and stars 
(scale-lengths hg = 25 pc and = lOpc, Mg — 10^ Mq and 

M, = 2.6 X IO^Mq, respectively; stellar disk with vertical sech^ 
profile and dispersions such that 2=1, gas disk initially thermally 
supported with hjR — 0.2). The initial surface densities of the gas 
and stellar disk are thus ~ 10^ M© pc“^ ~ 10 g cm“^. 
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Stellar Feedback, AGN Accretion, & Obscuration 3 



Figure 1. Morphology of the gas in a standard simulation, in face-on (x,y; left), side-on (x,z; middle), and cylindrical {R,r, right) projections. The time 
3Myr since the beginning of the simulation) is 150 (8) orbital periods at Ipc (10 pc). Brightness encodes projected gas density (increasing with density; 
logarithmically scaled with a 6dex stretch); color encodes gas temperature with blue material being T < lOOOK molecular gas, pink ~ 10^ — 10^ K warm 
ionized gas, and yellow > 10^ K hot gas. Top: Simulation with stellar, but no AGN feedback (no_BAL in Table [^. A multiphase disk forms; it is mostly 
molecular inside the central ~ 200pc, with heating by HII regions very localized to small ionized “bubbles” and heating by SNe restricted to low-density 
regions where it can vent vertically. The central ~ 10 pc develops a stellar+gas accretion disk dominated by m = 1 modes. Bottom: Same, with broad- 
absorption line winds (vSOOO). The winds blow out a polar cavity and generate an expanding shell in-plane, with occasional dense clumps sinking through to 
the center. Feedback eventually evacuates the entire nuclear region. 


This is chosen, based on the survey in [Hopkins & Quataert] 
( |2010a| ), to provide a “best case” for extremely rapid BH growth 
and quasar-level fueling. It is motivated by large-scale simulations 
of major mergers which produce dense, torus-like structures and 
high accretion rates ([Hopkins & Quataert|2010b||2011b|a[|Hopkins| 


[et al.[[2012a|b[ [Hop^ns[ 2012) , as well as at least some observa¬ 

tions indicating the presence of powerful AGN in nuclear starburst 
“cusps” even in galaxies which may not be experiencing extended 
star f ormation ([Diamond- Stanic & Rieke[[2012[ [Mushotzky et al.[ 


|M4l|Alatalo et al.|2015t.' 

The initial gas disk contains : 


; 0.6 X 10^ particles; the initial 


gas particle mass is 20Mq. We consider a limited resolution 
comparison in Appendix The force softening for the BH, gas, 
and star particles is set to e = 0.02 pc, with minimum SPH smooth¬ 
ing length = 0.1 times this. We note that all simulations employ 
the more sophisticated formulation of artificial viscosity described 
in [Morris & Monaghan[ (r997[ ), which greatly reduces numerical 
dissipation away from shocks relative to earlier implementations 
(see e.g. [Rosswog et al.|2000[[Price|2008[ l. 


2.2 Cooling, Star Formation, & Stellar Feedback 

Gas follows an atomic cooling curve with additional fine-structure 
cooling to 10 K. Metal-line cooling is followed species-by-species 
for 11 tracked species as in [Wiersma et al.[f2009a|bj t. The enrich¬ 
ment for each species is followed with the time dependent metal 


fiux directly attached to the mass, momentum and energy fiux from 
stellar winds and SNe Types la & II (see [Hopkins et al.[[2012c[ 
[2013b| l. 

Star formation is allowed only in dense, m olecular, self- 
gravitating regions above n > lO'^cm”^. We follow Krumholz & 
[Gnedin[ ( [2011j ) to calculate the molecular fraction in dense gas 
as a function of local column density and metallicity, and allow SF 
only from molecular gas. Following [Hopkins et al.[ ( [20T3c] ), we also 
restrict star formation to only gas which is locally self-gravitating, 
i.e. has a = 5v^{5r) (5r/Gmgas(< 5r) (1/4) || V(8)v|p/ (Gp) < 1 

(where the limit is taken as the “averaging radius” 6r vanishes, al¬ 
lowing a to be calculated in a resolution-independent manner only 
as a function of local properties). Gas which meets all of these cri¬ 
teria forms stars at a rate p* = pmoi/te (i.e. 100% efficiency per 
free-fall time). As shown in [Hopkins et al.[ ( |2013c[ ), the molec¬ 
ular criterion is not especially important in galaxy centers since 
most of the dense gas is molecular already, but the self-gravity 
criterion is important for small scales around black holes, where 
any simple constant-density threshold for star formation fails to ac¬ 
count for the radially-dependent tidal forces. Even in these regions 
however, the role of these criteria is primarily to determine where 
stars form; [Hopkins et al.[ ( |2013c[ ) and a number of other studies 
have shown that the total star formation rate, once fragmentation 
and stellar feedback are resolved, is set by stellar feedback, and 
is largely insensitive to details of both cooling and star formation 
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Figure 2. Face-on morphology of the gas in the additional simulations from 
Tablejl] as Fig.[T] at the same time (and same scale). Top: Compton heating 
(v5000_C) and isotropic vs. disk-planar winds (v5000_iso) have no visible 
effects (the differences are consistent with stochastic variations run-to-run). 
Middle: Lowering (v5000_loP) or raising (v5000_hiP) the momentum¬ 
loading of the winds leads to smaller/larger bubbles after the initial accre¬ 
tion event, although these in turn alter the subsequent accretion rate (see 
§ 1^. Bottom: Lowering the wind velocity (at fixed momentum-loading; 
v500) has no significant effect (although the shocked gas is colder, as ex¬ 
pected). Raising the wind velocity to 30,000kms“^ (v30000) creates the 
most hot gas (as expected); however, this gas appears to mostly vent out 
from the central regions, so a nuclear disk similar to the no-feedback case 
re-forms within ~ lOpc which drives a higher accretion rate at later times. 


prescriptions ( see|Saitoh et al. 12008}[Hopkins et aL|2011||2012d|c| 
|2013b|a||2015[|Agertz et al.|2013|l. 

Once stars form, feedback is included in the form of radiation 
pressure (UV, optical, and IR, allowing for multiple-scattering), 
stellar winds (fast, young star winds and slow AGB winds), SNe 
(types la and II), photo-ionization and photo-electric heating. Ev¬ 
ery star particle is treated as a single stellar population with an 
age based on its formation time and metallicity and mass inher¬ 
ited from its parent gas particle. Feedback includes the relevant 
mass, metal (with 11 separately tracked species), momentum, and 
energy injection to the neighboring gas; all of the relevant quanti¬ 
ties (stellar luminosities, spectral shapes, SNe rates, wind mechan¬ 
ical luminosities, yields) for the mechanisms above are tabulated 
as a function of time directly from the stellar population models 
in STARBURST99, assuming a |Kroupa| ( |2002| ) IMF. For every SNe 
event (or every timestep for winds and single-scattering photon mo¬ 
mentum), the relevant energy, momentum, mass, and metals are 


deposited into the nearest gas particles surrounding each star parti¬ 
cle; long-range photo-heating and radiation pressure are treated in 
a simplified manner assuming spherically symmetric photon prop¬ 
agation from each star particle as an independent source. See |Hop-| 
[kins et all| < |2011|[2012d| l for details. The end result of this stellar 
feedback is a multiphase ISM with a broad range of densities and 
temperatures. 

2.3 Black Hole Growth & Feedback 

The simulations all include super-massive BHs. The BH is much 
more massive than the stellar/gas particles, so we do not need to 
artificially “force” the BH particle to stay in the center of the poten¬ 
tial, but let it move freely. We cannot, however, directly resolve the 
viscous accretion disk of the BH on scales 0.1 pc. We therefore 
simply assume that the BH immediately accretes any gas particle 
gravitationally bound to it (relative velocity less than escape) with 
apocentric radius (calculated from the particle position and veloc¬ 
ity relative to the BH) < 2.8e (the minimum Keplerian distance). 
The rate of particle accretion is capped at the Eddington limit. 

The BH radiates at a luminosity L — erMsHC^ (cr = 0.1 is 
assumed)Q The explicit details of the BH feedback implementa¬ 
tion are given in Appendix [A| we briefly summarize them here. 
Since quasars are believed to have high-velocity, near-planar winds 
driven off the accretion disk (e.g., [Murray et al.|1995| l, we assume 
that a fraction of the photon momentum drives a wind launched at 
the resolution scale around the BH from accreted gas. Specifically 
a fraction of any gas accreted is blown out as a wind with veloc¬ 
ity Vwind, planar with the inflow (by launching particles directly at 
the accretion radius with this velocity). Two parameters define the 
wind, the mass-loading and velocity; this is equivalent to speci¬ 
fying the momentum-loading (/iwind = ^pL/c) and energy-loading 
(£^wind = tie L) of the wind. Values for the simulation parameters are 
in Tabled 

We also include Compton heating & cooling from the radi¬ 
ation field. Following [Sazonov et al.| ( |2004| ), this can be approx¬ 
imated with a nearly obscuration-independent Compton tempera¬ 
ture of Tcompton ~ 2 X 10^ K. Wc add the appropriate Compton rates 
to the standard cooling function (with a limiter following |Faucher-| 
[Giguere & Quataert|2012| to account for rate-limiting by Coulomb 
collisions at the high temperatures that can obtain in strong shocks). 

3 RESULTS WITH STELLAR FEEDBACK, BUT NO 
BLACK HOLE FEEDBACK 

Fig. {top row) shows the morphology of the high-resolution no 
AGN feedback run at a typical time after a few orbital periods, 
when the system has reached an approximate statistical steady state 
(Fig. [^compares the additional simulations in Table [3- Figure!^ 
shows the star formation rate as a function of time for simulations 
with and without AGN feedback {top panel) and two versions of 
the Kennicutt-Schmidt relation describing the star formation law 
for these nuclear-scale simulations {bottom two panels). Note that 
in the simulation with only stellar feedback, there is an initial burst 
of star formation but after a few Myr, the star formation rate settles 
into an approximate steady state at ~ 1M© yr“^ within ~ 1 kpc. 
The image in Fig. is shown in the latter phase. Within < 10pc, 
stellar feedback alone does clear most of the gas after a few Myr; 
this is recycled in a small-scale fountain on a similar timescale. 


^ We also describe in Appendix a model which imposes a spectrum of 
sub-grid time variability in the accretion rates; however this has no signifi¬ 
cant effects on the time-averaged results here. 
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lOg(Egas) [M®pc^] 



log( a ) [ Mgyr' kpc'^ ] 

Figure 3. Top Panel: Star formation rate within 10 and 500 pc regions for 
simulations with (vSOOO) and without (no_BAL) AGN feedback. Middle 
and Bottom Panels: Location of the same simulations on the Kennicutt- 
Schmidt relations at different times. The star formation rate surface density 
and gas surface densities are averages within 10 pc and the rotation rate in 
the bottom panel is also measured at 10 pc. The observations in the mid¬ 
dle panel (dashed line ±0.5 dex is from |Narayanan et^j2012^ ’s variable 
Xco model) are based on a range of galaxies, not just galactic nuclei, but 
nonetheless provide a useful point of comparison. The star formation effi¬ 
ciency per dynamical time evolves significantly with time during the sim¬ 
ulation, with a relatively high star formation efficiency in the burst of star 
formation at early times followed by a more prolonged period of lower star 
formation efficiency. 



0.1 1 10 100 



0.1 1 10 100 
R [pc] 


Figure 4. Time-averaged structural properties of the simulations. Top: 
m = I mode amplitude |flm=i I m the cold molecular gas, as a function of 
radius. With no BAL feedback the large spiral modes are visible here; with 
BAL winds the order-unity asymmetries introduced by the AGN wind im¬ 
pacting the ISM dominate. Bottom: Gaussian disk scale-height (h/R) versus 
radius. With no BAL winds, a modest h/R ~ 0.1 — 0.2 ~ \am=\ \ is sup¬ 
ported by the combination of stellar feedback and gravitational instabilities. 
With BAL winds, h/R is greatly enhanced because there is little gas and it 
is often dominated by escaping/venting polar winds. Even in the latter case, 
the scale-height of the cold rotating gas remains modest, similar to that in 
the non BAL wind simulation (see Fig.[^. 


The dynamics of these small-scale burst-quench cycles is explored 
in more detail in |Torrey et al.|p01^ . 


3.1 Black Hole Accretion 

Fig. □ shows that the gas disk exhibits strong non-linear m — I 
spiral wave and eccentric/lopsided disk modes, which are visible in 
spite of the inhomogeneous structure of the ISM. Using simulations 
on similar spatial scales but with a much less realistic model of the 
ISM, [Hopkins & Quataert| ( |2010a| ) showed that non-linear m = 1 
modes generated by stellar-gas interactions dominate the angular 
momentum transport in galactic nuclei at and inside the BH sphere 
of influence. However, that study was limited by the assumption 
that the ISM gas was described by a simple “effective equation of 
state” (i.e. no resolved phase structure, winds, or turbulence, simply 
single-phase gas with a barytropic non-thermal pressure following 
[Springel & Hernquist||2003| l. We confirm their result here with a 
much more realistic ISM model. 
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Figure 5. Top: BH accretion rate vs. time. To make systematic differences 
clear we smooth the rates in a ~ 2 x 10^ yr window (see Appendix[^for an 
un-smoothed case). AGN feedback suppresses the BH accretion rate rela¬ 
tive to simulations without AGN feedback, with higher momentum-loading 
in the input AGN wind (higher rjp) leading to lower Mbh- Middle: Total 
momentum flux in the galaxy-scale outflow at > lOpc in each model, vs. 
time. The models with BAL winds all equilibrate at broadly similar outflow 
momentum flux. This explains why higher-77p models adjust to have lower 
Mbh- Bottom: Corresponding mass-outflow rate of the wind at > 100pc. 
Simulations with AGN feedback have dramatically larger outflow rates 
from the nuclear region relative to the simulation without AGN feedback. 


Fig.|^(?o;? panel) plots the m = 1 mode amplitude^ versus 
radius for the simulations with and without AGN feedback. For 
the simulation without AGN feedback, the m—\ mode amplitude 


^ Mode amplitudes are measured in the gas surface density as 


\am{Rd)\ 


I S(/?, (j)) exp (imcj)) d(f)\ 

4')d4> 


( 1 ) 



Figure 6. Inflow/outflow rate vs. radius in the model with stellar feedback 
alone, averaged over time for a few dynamical times. Negative values (out¬ 
flow) are dotted, positive values (inflow) are solid (absolute value plotted 
for the sake of a logarithmic projection) “We compare several analytic ac¬ 
cretion rate models (§ [^. The “gravitational torques” estimator (resonant 
exchange between gas+stellar gravitational instabilities) is accurate within 
a factor ~ 3 at all radii < 100pc, and correctly predicts the major sign (in¬ 
flow/outflow) changes. Spherical accretion models fare poorly: the “pure 
Bondi” estimator gives ~ lO^MQyr”^ too large to fit on the plot; the 
“modified Bondi-Hoyle” estimator over-predicts by ~ 2 — 4dex. “Ballistic 
accretion” from turbulence fares poorly in the opposite manner (predicting 
^ 10“^MQyr“^ at R < 40pc). “Gravito-turbulent viscosity” is dimen¬ 
sionally reasonable, but under-estimates M by factors of ~ 5 — 50 near the 
BH radius of influence, where gravitational torques are most prominent, and 
does not capture the sign information. Around ~ 0.1 pc, resolution effects 
from our discrete particle number become important (some predictions drop 
precisely because the BH is accreting particles). 


found here ~ 0.1 is similar to that found in [Hopkins & Quataert| 
( |2010a| )’s simulations with gas fractions > 0.5. This suggests that 
the mode excitation and saturation physics is at least broadly sim¬ 
ilar in spite of the more dynamic multi-phase ISM present in our 
simulations. 

Fig. shows the black hole accretion rate|^ the outflow rate 
from the galactic nucleus, and the total momentum flux in the 
outflow as a function of time|^The simulations clear ly And large 


inflows up to ~ Mq yr“^ to the central < 0.05 pc. Hopkins & 
IQuataertj l |2011a| ) derive an analytic approximation for the inflow 
rate through each annulus for inflows driven by strong gravita¬ 
tional torques and resonant angular momentum exchange between 


^ To highlight the systematic differences between runs, we boxcar-smooth 
each curve in a time window of 2 x 10^ yr. In Appendix we show that 
there is variability on all resolved timescales in the simulations (as also 
seen by jNovak et al.|201 ijjGan et al.|2014[|Dubois et al.|2014) , and we even 
consider a model for un-resolved time variability. However, we caution that 
some of the resolved small-timescale variability is almost certainly artificial 
here (owing to the assumption that individual gas particles are accreted dis¬ 
cretely) - there are ~ 10^ — 10^ such accretion events per simulation. The 
BH accretion rate would be likely smoothed out if these were accreted into 
a viscous disk, which then accretes onto the BH. 

^ The net inflow rate at radius R is given by M = AR~^ f dMgas vr in an 
annulus. The outflow rate is the same integral, but only over dMgas where 
Vr > 0. The rates are time-averaged in each annulus (which also removes 
the spurious radial velocity contribution from e.g. stationary modes). Be¬ 
cause of finite bin-widths the inflow rate can change sign discretely from 
bin-to-bin. 
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gas and stars. For modes with complex potential ^a{R) and pattern 
speed cj = Qp + / 7 , this is: 


M=Egas/?"n|^ 


'a\\n 

A U 


+ d\nVc/d\nR\ 


( 2 ) 


with 5(ct;, $a) a phase function and F an order-unity amplitude 
correction derived in [Hopkins & Quataert| ( |2011a| ), which can be 
measured directly in the simulations. For an m = 1 mode in a quasi- 
Keplerian potential, this is approximately M ~ — |a| Egas R^ Ft (with 
\a\ the mode amplitude). 

Fig. E compares equation]^ to the simulation inflow/outflow 
rates; the agreement is reasonable, particularly given that the an¬ 
alytic result was derived under the assumption of smooth (non- 
turbulent) gas flows. Note that divergence between models at < 
0.1 pc is primarily a consequence of our discrete accretion model 
actually removing particles at this radius, and should be considered 
with caution. 

Fig.j^also compares the inflow rate in our simulations to four 
alternative proposed accretion rate estimators, none of which does 
as good a job of reproducing the simulation results. (1) Bondi: 
MBondi ~ dTrG^MgHPgascE^- TMs ovcr-prcdicts the accretion rate 
by an enormous factor ~ 10^ as most of the gas is cold and 
molecular, supported not by pressure but by angular momentum. 
(2) Modifled Bondi-Hoyle: Mmbh ~ 47rG^Menc(< + 

(7as -bh)) This allows for the fact that mass outside the BH 
itself (e.g. the bulge or nuclear cluster) should, when viewed from 
gas at large enough distance, act as a point mass in the same 
way; it also allows for super-sonic relative motion of the gas and 
BH. This is dimensionally closer to what we measure than [1], 
but given the low c^, it amounts to assuming all gas is in free- 
fall (neglects angular momentum) - i.e. it is quite similar to sim¬ 
ply taking Mmbh ~ Mgas/feee-faii - and it over-predicts M by fac¬ 
tors ~ 1000. Note that replacing Menc ^ Mbh, as in the standard 
Bondi-Hoyle formulation (used in Springel et al.||2005| Hopkins] 


|et al.|2006a||2005^|Di Matteo et ai.|2008 [Croft et al.|2()09| l only 

slightly decreases the discrepancy. (3) Ballistic Accretion: Mbaii ~ 
Itt Egas R^ ^ {Vc /cr) exp (—9 Vc^/16 cr^) (this corresponds to accre¬ 
tion of the randomly-populated low-angular momentum “tail” of 
highly turbulent flows from [Hobbs et al.|2011[ we generalize their 
formulae for accretion through each annulus). This disagrees with 
the simulations as well; dimensionally it gives M oc Mgas{R) F{R) 
but with a “reduction factor” {h/R)~^ Qxp{—0.56 [h/R]~^), 
which for /z//? ~ 0.1 — 0.3 found here is very small, so that there 
is very little ballistic accretion. (4) Gravito-turbulent viscosity: 
Mturb ~ 37racrgasEgasf^~^ whcrc a ~ 0.005 — 0.05 is the (cool¬ 
ing function-dependent) effective turbulent viscosity for a g = 1 
disk ([Gammie[[200r| [Thompson et al.[[?005[ [Debuhr et al.[[2010] 

201 1| ) ^This is dimensionally similar to the gravitational torques 

scaling but with free-fall slowed by a term a (h/R)^ instead of |a|; 
over some radii the two are comparable but the former decreases 
rapidly inside the BH radius of influence (implying accretion would 
be “throttled”) while \a\ can remain order-unity all the way to the 
true accretion disk (see[Tremaine[1995|[Bacon et al. [200 1| [Hopkins [ 
[& Quataert|20TT^[2010b[[Hopkins|2010jl~ 

The comparisons in this section are based on simulations with¬ 
out AGN feedback. In the presence of feedback, the net accretion 
rate onto the BH is determined by a competition between the in- 


^ In Fig. we adopt a = 0.05, close to the predictions from [Gammie[ 
( [2001) given the measured Mach number in the diffuse gas. Changing the 
value of a will systematically shift the normalization of the predicted curve. 


flow rate from large scales set by gravitational torques and the ef- 
flciency of AGN feedback at suppressing this inflow in the galac¬ 
tic nucleus. Our simulations explicitly resolve this competition and 
produce accretion rates a factor of ~ 10 lower than in simulations 
without AGN feedback (Fig.[^. For lower resolution galaxy-scale 
or cosmological simulations it is unclear what the best time aver¬ 
aged accretion rate estimator is to capture this competition between 
inflow by gravitational torques and AGN feedback; this merits fur¬ 
ther study in future work. 

3.2 Star Formation and Vertical Disk Structure 

Fig.[^(to/7 panel) shows the vertical scale height of the gas disk as a 
function of radius. The disk is in vertical equilibrium but the disper¬ 
sions are turbulent (much larger than thermal). As shown in Paper 
II on larger scales, stars form roughly until feedback can maintain 
Toomre 1 {h/R ^ Mgas(< /?)/Menc(< R)) and offset further 
collapse. At large radii this gives /z//?~0.2 — 0.3;atr~3 — lOpc 
this is h/R ~ 0.1 (dispersions ~ 20 — 70kms“^). The cylindrical 
image in Figure highlights the modest thickening of the disk at 
larger radii that is qualitatively analogous to that required in AGN 
"torus" obscuration models. As we describe in 0 this effect is 
much more dramatic in simulations with AGN feedback because 
feedback efficiently evacuates the polar region of gas. 

The bottom panels of Figure shows our simulations in two 
common versions of the Kennicutt-Schmidt relation. The star for¬ 
mation rate surface density and gas surface densities are averages 
within 10 pc and the rotation rate in the bottom panel is also mea¬ 
sured at 10 pc. The observations in Figure are best fits from 
[Narayanan et al.[ ^20l2\ based on a variable Xco factor. They are 
shown to provide a point of comparison, but include a range of 
galaxies, not just galactic nuclei. The time averaged star formation 
efficiency in Figure |^is broadly consistent with observations. The 
efficiencies evolve signiflcantly with time, however, with a rela¬ 
tively high star formation efficiency in the burst of star formation 
at early times followed by a more prolonged period of lower star 
formation efficiency. Perhaps most striking is that the star forma¬ 
tion efficiency per dynamical time decreases by nearly a factor of 
~ 10 during the course of the simulations without AGN feedback. 
Thus the decline in the star formation rate is not simply due to gas 
depletion but is also due to the decreasing star formation efficiency. 
Note that the duration of the simulation is comparable to the life¬ 
time of massive stars. Thus the stellar feedback that is effective for 
most of the duration of the simulation is that due to stellar radiation 
and stellar winds, since supernovae only start after ~ 3 Myr and 
have not had signiflcant time to operate. In addition, because the 
local dynamical time is short compared to the lifetimes of massive 
stars, the efficiency of stellar feedback depends primarily on the 
surface density of young stars, rather than the star formation rate. 
We explore the consequences of this for the “burstiness” of nuclear 
star formation and origins of the nuclear-scale Kennicutt-Schmidt 
relation in a companion paper ( [Torrey et al.|2016[ ). 

4 RESULTS WITH BLACK HOLE FEEDBACK 

We now consider the results of simulations with AGN feedback, 
focusing on the flducial v5000 run in which the AGN wind at small 
radii is injected with p — Ljc and F — 0.008 L. However, we con¬ 
sider variations in both the mass and momentum-loading as well, 
as outlined in Table [T] 

V\g.^{bottom) shows the gas morphology at a few Myr in our 
v5000 simulation; there is a clear dramatic impact of feedback on 
the gas, with the central ~ 30 pc relatively evacuated of gas by the 3 
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Figure 7. Column density distribution on sightlines to the BH in each sim¬ 
ulation. We integrate the column along 1000 sightlines to each BH at each 
time following [Hopkins et al.H2005dt , uniformly sampling the sky in solid 
angle, and show the distribution over all sightlines and times in each sim¬ 
ulation. Stellar feedback alone produces a relatively narrow range of very 
large columns. Simulations with BAL winds have evacuated polar regions 
with column densities < 10^^ cm“^ and overall broader obscuring column 
density distributions. The “clumpy torus” thus naturally arises from AGN 
feedback interacting with the large-scale ISM at ~ 10pc. 
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Figure 8. Median column density as a function of polar angle 0 (averaged 
over azimuthal angle cj)), calculated over all times as Fig. [7] for the simu¬ 
lations in Table [2 The disk mid-plane {cos 6 = 0) features Compton-thick 
columns, as expected, which decline towards the poles (cos^ = ±1). By 
clearing the central regions, feedback suppresses the column densities at 
all 6, but proportionally more towards the poles (where the smaller initial 
column makes complete evacuation easier) - this is needed to create opti¬ 
cally thin sightlines. In all the simulations, the 1 a scatter in column density 
at a fixed 0 is ~ 0.5 dex, owing to sub-structure in the gas (Fig.[^. Thus 
the ~ 95% inclusion range for v5000 (shaded) is larger than the system¬ 
atic offset between it and any other simulation with feedback (but not the 
no-feedback case). 


Myr time of these images. Qualitatively, our other simulations with 
feedback (Fig. resemble this case, but with subtle differences 
discussed below. 

The v5000 outflows are launched in the dense disk mid-plane. 
This drives an expanding shell in the disk plane, with gas piled up in 
a narrow ring/shell at the outer (radiative) shock where the winds 
are encountering the ISM. This is similar to |Faucher-Giguere &| 
IQuataert] pOI2| )’s models for galaxy-scale winds driven by AGN, 


though it is not clear if those models quantitatively apply because 
the hot shocked gas created by the AGN wind is not well-conflned - 
this may be a limitation of the small scales we are simulating (there 
is no full galaxy and halo into which the winds can propagate in 
these simulations). Indeed, out of the midplane, the entrained mass 
is modest so outflows coast or are accelerated by hot gas pressure 
Ailing the growing central cavity in the disk. 

The large impact of the AGN wind on the ambient gas has 
three closely related effects. First, it strongly suppresses the star 
formation in the galactic nucleus, by a factor of ~ 10 — 30 (Fig. 
1^. Secondly, it increases the net outflow rate from the galactic nu¬ 
cleus by a factor of ~ 10 — 30, to ~ M© yr“^ (Fig.[^. Finally, on 
longer timescales the BH feedback roughly regulates the BH ac¬ 
cretion rate. Speciflcally, the feedback momentum flux scales as 
p — r\pLlc\ balancing infall with feedback therefore implies a crit¬ 
ical value of /), so in equilibrium (L) . This scaling provides 

a reasonable approximation over a sufficient time average (Fig.[^, 
but the evacuation of the central regions clearly leads to very large- 
amplitude variability on ~ 10^“^ yr timescales. 

It is useful to directly compare the momentum flux in the 
galaxy scale winds (Fig.[^ with those injected at small radii in the 
AGN wind. These need not be the same if AGN feedback produces 
a bubble of hot gas that does work on the surrounding material, 
increasing the momentum flux in the wind (e.g., |Faucher-Giguere| 
|& Quata^|2012| t. For low input wind velocities the momentum 
fluxes in Fig. are comparable to that injected in the AGN wind 
at small radii, while for higher input wind velocities (in particu¬ 
lar, the v30000 simulation), there is a factor of few boost in the 
AGN wind momentum flux - although we inject a momentum flux 
~ L/c, the outflow momentum fluxes reach ~ lOL/c, comparable 


~ 1000kms“‘ ( 

|Borguet et al.| 

|2012 Cimatti et al.| 

|2013 Cicone 

et al.[[20141 [Harrison et al.[[2014|[Zakamska & Greene[[2014|[Za- 


shocked heated by the AGN wind is able to escape relatively easily 
along the polar direction. In a more self-consistent calculation, it 
is possible that the existence of large warps between the disk axis 
on small scales and that on large scales (e.g., [Hopkins et al.|2012b| l 
might act to better conflne the outflow. 

Note that the outflow rates of ~ 1 — IOMq yr“^ that we And 
are relatively modest, compared to many observations at >kpc 
scales (see references above). This, of course, owes in part to the 
limited region we are simulating (< 100 pc) - there simply is not 
much material in this volume for the winds to “sweep up.” On 
these scales directly, there are actually relatively few constraints 
on AGN outflow rates and velocities (most of the constraints above 
apply either to AGN accretion disk scales, or spatially-resolved out¬ 
flows on ~kpc scales). However, if we assume the outflow remains 
momentum-conserving, then for typical galaxy mass profiles we 
would easily expect it to entrain an order or magnitude or more 
mass as it propagates from <0.1 kpc to a few kpc. 

Figs. [7| 8 [ plots the column density distribution and dependence 
on viewing angle (both averaged over time and in various time in¬ 
tervals) for each of our simulations, both with and without AGN 
feedback. The model without AGN feedback predicts virtually no 
systems with column densities below 10^^ cm“^ even along polar 
sightlines, in stark contrast to observations. Even though the ISM 
is highly inhomogeneous on larger scales, a small, dense thick-disk 
or “halo” component surrounding the BH in the central ~ 0.1 pc is 
sufficient to produce these extremely high column densities even in 
the polar direction. The BAL winds have, however, an enormous 
impact on the column density distribution. This is not surprising 
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given their impact on the nuclear gas morphology. The polar re¬ 
gions are completely evacuated, giving a large fraction of sightlines 
that are fully un-obscured. The remaining sightlines follow a broad 
column density distribution, driven in part by the fragmentation and 
asymmetries seen in the expanding equatorial shells. The evacua¬ 
tion of the central regions out to some radius where h/R ^0.1 —0.3 
gives a canonical “torus-like” global morphology. This is particu¬ 
larly clear in the cylindrical image shown in the lower right panel 
of Figure 

4.1 Dependence on the Strength and Form of AGN Feedback 

In Table we outline a series of runs changing the energy and 
momentum loading of the AGN-driven winds. Here we compare 
their properties. 

First, we add Compton heating/cooling to our “standard” case 
(run v5000_C). Consistent with previous studies {|Ostriker et al. 
20T0||Ciotti et al.|2010|[Choi et al.|2012| 1^4] [2^1 [Park et al. 

2014| ) and observational constraints ( [Stern et 'ni||2015[ jChatter- 
jee et al.||2015| ), this has little or no effect on the properties we 
measure (see Figs. |2|7j ). This should not be surprising: all but 
the lowest-density gas in the simulations has cooling times much 
shorter than the Compton heating time. To the extent that Compton 
heating/cooling is important, it has been speculated that it may be 
important for cooling in the reverse shock of the BAL winds; how¬ 
ever, in the cases we simulate this is not expected to dominate even 
for a spherical blastwave (see jFaucher-Giguere & Quataert|2012] |, 
and “venting” rather than cooling dominates the escape of energy 
in the hot gas. 

We also consider a case (v5000_iso) where the winds from 
the BH are directed isotropically from the BH, as opposed to in 
the accretion plane. This is discussed in §|^ Figs. |2|^ demonstrate 
this produces only very small changes (within the range produced 
by purely stochastic effects) relative to our standard v5000. This 
is also consistent with previous studies ( jPebuhr et al.|2012] ), and 
follows simply from the fact that the distribution of accretion di¬ 
rections/angles (the fractional scale-height of accreting gas) is rel¬ 
atively large. 

Our v5000_hiP and v5000_loP runs keep the wind velocity 
fixed but increase/decrease the mass (and momentum) loading by 
an order of magnitude, respectively, relative to the BH accretion 
rate. Not surprisingly, a lower (higher) mass-loading produces ini¬ 
tially weaker (stronger) outfiows: the “bubble” in Fig.j^is smaller 
(larger), and this produces slightly higher (lower) column density 
sightlines in Fig.|^ This has relatively weak effects on the gas prop¬ 
erties outside the region being evacuated (Fig. |^. The SFR and 
Kennicutt-Schmidt relation are essentially indistinguishable from 
those in Fig.j^for our v5000 run, except higher/lower mass-loading 
translates to earlier/later suppression of the SFR in the immediate 
vicinity of the BH. The initially stronger outfiow (for the same ac¬ 
cretion rate) in v5000_hiP suppresses the BH accretion rate sub¬ 
stantially in Fig.|^ while the weaker outfiow in v5000_loP allows 
more rapid growth of the accretion rate compared to vSOOO. Inter¬ 
estingly, this ends up producing a very similar total momentum flux 
and outfiow rate in the winds at > 4Myr. This suggests that indeed 
we are seeing self-regulation when the BH injects sufficient mo¬ 
mentum into the medium to drive outflows that clear its vicinity. 
This critical asymptotic value is reached eventually in the differ¬ 
ent runs, and appears to only weakly depend on the feedback pa¬ 
rameters. However, how quickly the critical point is reached, and 
(correspondingly) how much the BH is able to grow, is strongly 
dependent on the initial feedback mass-loading. 

In our v500 and v30000 runs, we keep the momentum-loading 


fixed but vary the initial wind velocity from 500 — 3 x lO^krns-^ 
Decreasing the velocity to 500kms“^ leads to slightly higher av¬ 
erage accretion rates in Fig.[^ but the effect on the circum-nuclear 
structure (Figs. |2|4j ) and column density distribution (Fig. 1^, and 
outfiow rates (Fig. are weak. This is not surprising since in 
this limit, the winds are primarily momentum-conserving and so 
the absolute wind speed (at fixed momentum) does not qualita¬ 
tively change the dynamics. With ~ 3 x 10^kms“^ winds, however, 
we see interesting differences. The faster velocity produces more 
shock-heated high-temperature gas (unsurprisingly); the lower ini¬ 
tial wind mass-loading, however, also appears to lead to more effi¬ 
cient “venting” of the wind. This means there is less of a circum- 
nuclear “bubble” carved out in the cold/neutral gas in Fig. but 
rather more pronounced hot gas channels escaping. This in turn 
allows a thin nuclear disk to re-form (visible morphologically in 
Fig. 2 but also manifest in higher typical column densities in 
Fig. 71, which then produces higher BH accretion rates seen in 
Fig. 5 These higher infiow rates rates lead to a larger net wind 
momentum injection and outfiow rate. Some of these differences 
are similar to the behavior seen in the “thermal energy deposition” 
models of jChoi et al.| ( |2014l|2015| ); however, those simulations did 
not include the physics driving the multiphase structure of the ISM, 
and so could not capture the full magnitude of “venting” effects we 
see here. 


5 DISCUSSION 

We have used simulations with <0.1 pc resolution to study BH 
accretion and feedback in gas-rich nuclear disks around massive 
BHs accreting at quasar-like luminosities. Our calculations include 
an explicit treatment of star formation and stellar feedback, which 
produce a self-consistently inhomogeneous ISM. We model AGN 
feedback via Compton heating/cooling and high-speed accretion 
disk winds injected at small radii. 


5.1 The Role of Stellar Feedback 

Absent AGN feedback the properties of the gas disk inside ~ 100 pc 
are as follows. Gas cools efficiently and collapses in a mini- 
starburst until sufficient young stars are formed to maintain 2 ~ 1 
(mostly via radiation pressure-driven turbulence), leading to dis¬ 
persions ~ 20 — 100 kms“^ in a cold nuclear molecular disk. As in 
previous simulations which adopt highly simplified sub-grid mod¬ 
els of the ISM ( [Hopkins & Quata^|2010aj l, the disk develops 
large-amplitude m = 1 modes in gas and stars, and resonant angu¬ 
lar momentum transfer between gas and stellar disks drives rapid 
infiow of gas, with accretion rates of ~ 0.1 — 1 M© yr“^ at < 0.1 
pc. This agrees well with the analytic ( [Hopkins & Quataert|2011a[ l 
predictions for “gravitational torque”-driven accretion. In contrast, 
the Bondi-Hoyle, viscous, or ballistic accretion rate estimators fail 
to capture the simulation results and are not appropriate for the 
regimes simulated here, in which much of the gas resides in a rota- 
tionally supported (albeit geometrically quite thick) disk (Fig.[^. 

Stellar feedback does operate somewhat differently in galactic 
nuclei, as opposed to larger galactic radii, because the local dynam¬ 
ical time Q“Ms <Myr. Young massive stars are sheared into an 
un-clustered mass distribution (e.g. executing hundreds of orbits at 
~ 1 pc) before they explode. Rather than local, Jeans-scale clouds 
evolving independently, we should think of the disk as a coherently 
evolving, disky “star cluster” (see e.g. [Thompson et al.|2005| l. On 
longer timescales, this leads to episodic “burst-quench” cycles on 
small scales, studied in[Torrey et al.[([2016]). 
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5.2 The Role of AGN Feedback 

Our calculations demonstrate that high-velocity winds from the 
central < 0.1 pe with momentum fluxes ~0.1 — IL/c suggested by 
observations (e.g.,|Cimatti et al.|2013[|Tombesi et al.|2015[|Harri-| 

[son et al.|2014[|Zakamska et al.|2015| l have a dramatic effect on the 
circum-BH ISM. In particular, such winds can evacuate gas from 
the circum-BH disk (see Figs. [IE- This suppresses the star forma¬ 
tion rate and black hole accretion rate in the galactic nucleus by a 
factor of ~ 10 and enhances the gas outflow rate at ~ 100 pc by 
a comparable factor (Figs.[^&[^, also similar to observations in 
at least some systems with powerful on-going outflows ( |Shimizu| 
|et al.|2015l|Guillard et al.|20T5l|Alatalo et al.|2015| l. As expected, 
the amount of BH growth required to produce this level of feedback 
and evacuate gas from the central regions depends inversely on the 
mass and momentum-loading of the BH accretion-disk winds. The 
amount of hot gas generated and its “venting” depend on the initial 
wind velocities. Our simulations thus provide support for models 
in which luminous AGN signiflcantly disrupt the ISM of their host 
galaxies, at least on scales <100 pc. Our simulations also specif¬ 
ically support the hypothesis that luminous AGN may play a key 
role in driving galaxy-scale outflows from gas-rich galactic nuclei. 

In the plane of the circum-BH disk, the AGN winds deceler¬ 
ate as material is entrained into expanding rings/shells. In the po¬ 
lar direction, however, the galaxy-scale outflows powered by the 
AGN retain high velocities (~ 5 — 30 x lO^kms”^) as they reach 
kpc scales; although not isotropic, the opening angle for the high- 
velocity outflow is large (> 2/3 of the sky). A detailed compari¬ 
son with observations is outside the scope of this work, because we 
consider only one initial condition, and do not model the large-scale 
galaxy beyond 100 pc on which many AGN-driven outflows are 
observed (although a follow-up study designed to compare in detail 
with these observations is in progress). However, the broad range 
of velocities present simultaneously within the same system is eon- 
sistent with outflows observed outside of accretion disk scales in 
molecular and ionized gas (see e.g.|Tom besi et al.|2015[ |Harrison| 

|et al.|2014||Zakamska & Greene|2014[|Alatalo|2015| l. Siinilarly, we 

identify outflowing material across a broad range of temperatures 
and spatial scales, from molecular to > 10^ K gas, and from scales 
of ~ 0.1 — 1000pc. Much of the material in the disk plane is ac¬ 
celerated to low velocities and will not escape, but entrains a large 
mass (comparable to or larger than starburst-driven winds). At least 
some energy of shocked AGN winds is converted into work, allow¬ 
ing the outflow momentum to reach ~ lOL/c in some cases (again, 
qualitatively consistent with observat ions; see|B orguet et al.|2012[ 

Cimatti et al.|2013[|Cicone et al.|2014[|Harrison et ^.|2014[|Zakam-| 

ska & Greene|2014| ). This is dietated largely by the geometry of the 
surrounding ISM, rather than the AGN wind at small radii. In par¬ 
ticular, simulations with isotropically directed AGN winds on small 
scales give similar results to our default calculations that utilize pri¬ 
marily planar winds (this highlights that once the AGN wind shocks 
the gas follows the “path of least resistance” in the polar direction 
independent of exactly how the wind is initially directed). Under¬ 
standing whether the outflows we And will be conflned or halted by 
the galactic ISM or will continue to escape out of the galaxy will 
require galaxy-scale simulations. It is important to stress that the 
present caleulations are not well-suited for addressing this question 
because our idealized initial eonditions do not have, e.g., a gaseous 
halo or the nuclear warps/mis-alignments seen in both simulations 
and observations of galactic nuclei. 

In our calculations, AGN-driven outflows also have a dramatic 
impact on obscuration of the AGN itself. AGN winds evacuate the 


polar region to allow a fully un-obscured view of the BH. AGN 
winds thus self-consistently produce a torus-like morphology (see, 
in particular, the lower right panel of Fig. [^. Quantitatively, we And 
a broad column density distribution from 10^2 _ 1026 cm■^ in 


reasonable agreement with observations (e.g. 

Malizia et al.|2009| 

ITreister et al.||20091 |Risaliti et al.||19991 |Bur 

Ion et al.||2011| and 


references therein). The inhomogeneous nature of the ISM also in¬ 
evitably introduces large (~ 1 dex) variation in obscuring columns 
even at roughly flxed polar angle - similar to observational sug¬ 
gestions of “elumpy” torii {[Risaliti et al.|2002[[Mason et al.|2006 
Sanchez et al.||2006[ [Nenkova et al.||2008[ [Ramos Almeida et al.| 

2009|[H5nig & Kishimoto|20T0||Deo et al.|2011| |. 

We And that Compton heating/cooling from the AGN pro¬ 
duces weak effects on these scales, consistent with previous studies 
dChoi et al.|20T2l[M4l[M^|Park et al.|2014| ). 

5.3 Future Work: Other Scales & Forms of Feedback 


These simulations are a first exploration of the interaction between 
AGN and stellar feedback on scales between the BH accretion disk 
and galaxy. We focused on these scales beeause they are relatively 
unexplored and yet critical for understanding BH growth and the 
impact of AGN winds and radiation on the ISM. It is, however, also 
clearly important to extend our models to cover a broader range of 
spatial scales. On smaller scales, understanding the origin of AGN 
winds and radiation and their “escape” from the accretion disk is 
critical for setting the magnitude and geometry of AGN feedback 
on pc scales. On galactic scales, we need to understand how the out¬ 
flows found here interact with the galaxy ISM on long timescales: 
in particular how this changes galaxy star formation histories and 
regulates future episodes of BH inflow. This is necessary to de¬ 
termine the effects of feedback on BH-host galaxy correlations. 
Also, many observations of galaxy-scale, AGN-driven winds sug¬ 
gest large momentum-loading in the winds, with up to ~ lOL/c 


(see Sturm et al.|2011[[Faucher-Giguere & Quataert[2012[|Cicone| 

|et al.|2014| l. It will be particularly interesting to see whether these 

observations are consistent with a model in which this large mo¬ 
mentum flux is generated on accretion disk scales (our v5000_hiP 
model), e.g., by super-Eddington accretion, or whether they require 
additional large-scale effects such as confinement (and buildup of 
a pressure-driven bubble) of radiation or hot shocked gas in the 
galaxies ISM. This could, e.g., be produced by misalignment be¬ 
tween the nuclear scale disk and the galaxy ISM as a whole. 

The two AGN feedback mechanisms we have studied here 
(fast AGN winds and Compton heating), are by no means exhaus¬ 
tive. In future work, we will extend this to include radiation pres¬ 
sure on both narrow lines and dust, photo-heating, and the effects 
of relativistic jets, all of which can act directly on gas both on 
the scales we model here but also on much larger scales up to 
> 100kpc. We will also study the effects of different initial con¬ 
ditions (e.g. gas fraction and disk-to-BH mass ratio). 
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APPENDIX A: BLACK HOLE FEEDBACK 
IMPLEMENTATION 

Al Broad Absorption Line Quasar Winds 

Bright quasars often have BAL winds with velocities of ~ 1000 — 
30000kms“^ We model these in the most direct manner possible: 
when a gas particle is accreted, a fraction /acc is actually assumed 
to accrete on the BH while the remaining 1 — /acc is blown out as a 
BAL wind with velocity vbal- There is both observational ( |Schmidt| 
|& Hines|1999[|Ogle et al.|199^ and theoretical ( [Murray et al.|1995| ) 
evidence that BAL winds are approximately planar (in or slightly 
out of the accretion disk plane). Assuming that the angular mo¬ 
mentum vector of the small-scale accretion disk is correlated with 
that of the sub-pc accreting material, then this corresponds to di¬ 
recting vbal along the radial vector R = ri — Tbh from the BH. On 
an accretion event we therefore take (for the accreted gas particle) 
Mi ^ (1 — /acc) m, apply the “kick” v/ ^ V/ + vbalR, and hold m 
(internal energy per unit mass) constant. In model v5000_iso we in¬ 
stead assign the wind direction randomly. Previous work has shown 
that randomly directed winds yield results similar to planar winds, 
but require somewhat larger wind momentum fluxes to achieve the 
same feedback on the ambient gas ( jPebuhr et al.|2012] l. 

Two parameters must be chosen: the initial outflow mass load¬ 
ing 13 = Mbal/Mbh = (1 - face)/face and velocity vbal; obser¬ 
vations and theoretical mo dels suggest values of order ^ 1. 
Vbal ~ 10"^kms ^ (see e.g. iMoe ^ al!]|2009[ [Punn et al.||2010| 


[Hamann et al.|20lf||Borguet et al.|2012| and references therein) 

Equivalently, we can translate these parameters into the wind 
momentum and energy-loading. Since BAL winds are believed to 
be driven by line radiation pressure in the accretion disk, the avail¬ 
able momentum flux is p — L/c, where L = CrMBuc^ is the lumi¬ 
nosity (with Cr ~ 0.1 the radiative efficiency). The “initial” wind 
momentum and energy are MbalVbal and 0.5Mbalv|al respec¬ 
tively; thus the energy and momentum-loading are 




Vbal 


30,000kms“ 

Vbal 

30,000kms‘ 


r) 


t)' 


(Al) 

(A2) 


Note for ? 7 p = /3 = 1, we recover the canonical ije ~ 0.05 adopted 
in previous simulations with purely thermal AGN feedback (e.g.[Pi| 
[Matteo et al. [2005[[Hopkins et al.|2005^ . 

A2 Compton Heating/Cooling 

The radiation held of the BH will also Compton heat/cool gas in 
its vicinity. As discussed in [Sazonov et al.[ ( [2004[ [2005[ ), this ef¬ 
fect is nearly independent of obscuration: Compton heating is en¬ 
tirely dominated by photons with energies ::$> lOkeV (for which we 
can usually safely ignore obscuration) and Compton cooling by the 
bolometric luminosity in lower-energy photons (re-distributed, but 
not, in integral, altered by obscuration). As such even Compton- 
thick columns result in factor < 2 changes in the heating/cooling 
rates. We therefore neglect obscuration and assume the radiation 
held is isotropic, so that the X-ray/bolometric flux from the AGN 
on all particles is given by Fx = Lx/^ ttj^, with Compt on tempera¬ 
ture ^2x 10^ K as calculated in 


Sazonov et al. 


( 2004\ for a broad 


range of observed QSO SEP shapes]^ In the cooling function. 


^ We propagate this flux through the gravity tree, since it follows an 
inverse-square law when we can neglect obscuration. This makes it triv- 


© 0000 RAS, MNRAS 000, 000-000 





























14 Hopkins et al. 



Time [Myr] 

Figure Bl. BH accretion rate vs. time as Fig. for our “no BAL” sim¬ 
ulation. We show the rate smoothed in a rolling 2 x 10^ yr window as in 
Fig. [ 5 ] (thick dashed), and smoothed in a 10^ yr window (solid). There is 
variability on all resolved timescales; however, some of this owes to numer¬ 
ical effects, namely the assumption of instantaneous accretion of discrete 
gas particles. In Appendix]^ we describe a model allowing variability even 
on unresolved (arbitrarily small) timescales. 


add the appropriate Compton heating and cooling terms [^Although 
Compton cooling depends explicitly on the free electron fraction, 
for the photon energies dominating heating (much greater than the 
ionization energy of hydrogen), we can safely approximate Comp¬ 
ton heating of bound electrons as identical to free electrons (see 
e.g. [Basko et al.|1974||Sunyaev & Churazov|1996| >. 

Finally, as shown in Faucher-Giguere & Quatae^ ( |2012j ), 
some care is needed at the highest temperatures: if the timescale 
for Coulomb collisions to transfer energy from ions to electrons 
is longer than the Compton or free-free cooling time of the elec¬ 
trons, this is the rate-limiting process and a two-temperature plasma 
develops. We therefore do not allow the Compton-i-free-free cool¬ 
ing rate to exceed the Coulomb energy transfer rate between ions 
and electrons calculated for an ion temperature T in the limit 
where the electrons are efficiently cooling Te T (see jSpitzerj 
|1962[[Narayan & Yi|1995j ). It is important to note that AGN wind- 
shocked electrons are generally non-relativistic: either immediately 
post-shock (where most energy is in protons, with electron tem¬ 
perature Te ~ Tp(me/mp) ~ 1.3 x 10^ K(vshock/30,000kms“^)^), 
or in later stages when competition between Compton cooling and 
Coulomb heating regulates the temperature. 


APPENDIX B: VARIABILITY ON UN RESOLVED 
TIMESCALES 

We are able in these simulations to follow infiows to sub-pc scales. 
This, coupled to the fact that our accretion model is by defini¬ 
tion discrete (i.e. particles are swallowed individually and instan¬ 
taneously accreted onto the BH) leads to large variability on very 
small timescales, illustrated in Fig. jBlj However, there are still sev¬ 
eral orders of magnitude between these scales and the BH event 


ial to apply the appropriate flux to arbitrary particle numbers, geometries, 
and numbers of black holes. 

^ As is standard, cooling is solved implicitly within this function in the 
regime where the heating/cooling times are short compared to the particle 
timesteps. 
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Figure Cl. Disk structure as Fig. but comparing different numerical 
methods (s ee §[C| our “pressure-entropy” formulation of SPH developed 
in |Hopkin's|(2013t , which performs well in tests of fluid mixing instabilities 
while maintaining good conservation; a “density-entropy” formation (the 
“standard” GADGET SPH method); a run with a simplifled artifleial viscos¬ 
ity prescription and SPH smoothing kernel; and a run with the “pressure- 
entropy” formulation but ten times as many particles. The results appear 
robust to these variations. 


horizon, spanned by the Shakura-Sunyaev accretion (a) disk. Em¬ 
pirically, AGN exhibit variability on all observed timescales, cor¬ 
responding to these unresolved spatial scales. Although we cannot 
resolve these scales, we can make a crude estimate of the effects 
of this variability by including a sub-grid power-spectrum of lu¬ 
minosity fluctuations and integrating over this to obtain the (mod¬ 
ified) momentum flux in every resolved simulation timestep. We 
quantitatively implement this following the prescription in jHopkinT] 
|& Quata^ ( |20 11 af , integrating over a power spectrum with equal 
logarithmic power per logarithmic time interval, from the minimum 
resolved timestep down to the orbital time at the innermost stable 
circular orbit for a non-rotating BH. Performing such an experi¬ 
ment, we find almost no effect on our conclusions. Given the re¬ 
solved dynamic range in the simulations, this additional variability 
occurs on extremely small timescales compared to the dynamical 
times of the outflow - the timescale over which feedback deter¬ 
mines the equilibrium accretion rate. As such, other than adding 
the chosen random variance to the lightcurve on small timescales, 
this introduces (relatively) little dynamical effect. 

APPENDIX C: NUMERICAL TESTS 

We now consider some tests of the robustness of the numerical 
methods used here. Figs. |C1|C^ repeat Figs. |4|5| but with varied 
numerical prescriptions. Our default simulations use the “pressure- 
entropy” SPH formulation described in [Hopkin^ pO 13 j ) , which is 
shown there to give dramatically improved results on in situations 
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to the asymptotic level of outflow occurs more rapidly at high resolution 
and more slowly with the modifled viscosity and kernel, but this is a result 
of the non-equilibrium initial conditions. The variations in the late-time ac¬ 
cretion rates owe to the chaotic bursts as individual cold gas clumps sink to 
the center, so we only expect statistical convergence in this stage. 


with fluid mixing around contact discontinuities (e.g. the Kelvin- 
Helmholtz and Rayleigh-Taylor instabilities) while retaining excel¬ 
lent conservation properties, and includes a number of additional 
improvements to the treatment of artiflcial viscosity (see |Cullen & 
Dehnen|r2010| ), SPH smoothing kernel accuracy ( [Dehnen & Aly 
2012| ), and timestep communication relevant for treating extremely 
high Mach-number shocks ( [Saitoh & Makino|2009[|Durier & Dalla| 
|Vecchia|2012| |. 

To test whether these subtleties may be strongly influencing 
our results, we re-run our standard vSOOO simulation instead using 
a “density-entropy” SPH formulation, as in jSpringel & Hernquistj 
( |2002j ) (the “standard” GADGET formulation of the SPH equations- 


of-motion). 0This produces a “surface tension” term at contact dis¬ 
continuities that suppresses some fluid mixing instabilities, which 
has been the subject of much discussion in the literature (see jAgert^ 
jet al.|2007l|Read & Hayfleld|2012| and references therein). We also 
re-run the simulations with the pressure-entropy formulation, but 
adopting the much simpler and more numerically dissipative con¬ 
stant form of artiflcial viscosity from jGingold & Monaghan| ( [l983j ) 
(which can signiflcantly alter the behavior in sub-sonic turbulence; 
see |Price||20f^ , and a greatly reduced-accuracy SPH smoothing 
kernel (a 32-neighbor cubic spline, as opposed to our standard 
128-neighbor quintic spline). Together these variations produce the 
range of numerical effects which span the major SPH-grid code 
differences often discussed in the literature (see references above 
andjPrice & Federrath|2010[|Bauer & Springel|2012[|Sijacki et al.j 

|2012| l. We also run a standard resolution test, increasing the number 

of particles by a factor of 10. 


We see very little difference in the results in Figs. |CHC2| 
Likewise there is relatively little difference in the column density 
distributions and gas morphologies. There are some slight differ¬ 
ences in the phase diagrams, which correspond to the degree of 
fluid mixing along phase boundaries (the quantity most affected by 
these differences), but it is mostly at low temperatures where it does 
not have signiflcant dynamical effects. This probably relates to the 
fact that in the numerical comparison studies discussed above, it 
is generally well-established that different methods agree well in 
the regimes of super-sonic turbulence and/or systems with domi¬ 
nant external forces. Moreover all of these changes preserve good 
energy and linear and angular momentum conservation, so to the 
extent that the outflow is primarily a simple momentum or energy- 
conserving “piston,” and the steady-state stellar feedback is the re¬ 
sult of momentum input balancing runaway collapse (see jHopkinsj 
jet aL|20lT]|2012d| ), our conclusions should be robust. 


In our resolution tests, we see quite good agreement in the BH 
accretion rate and wind momentum, up to the level of stochastic ef¬ 
fects (random differences between simulations). In fact re-running 
our standard model with different random number seeds for the 
placement of the initial particles in the disk leads to comparable 
variations. In the outflow rates, we see the largest differences be¬ 
tween runs, at times < 3 Myr. These also vary the most dramati¬ 
cally owing to purely random effects; at these times (order one dy¬ 
namical time at ~ 100pc), the outflow rate measured at 100pc has 
just begun to develop, so even small differences in the dynamics 
can appear signiflcant. The high-resolution case develops a wind 
more quickly because fragmentation and star formation in the disk 
at these radii, being better resolved, proceeds faster, making the 
disk thicker and putting more material into a stellar-feedback driven 
fountain or wind which the BH wind is then able to entrain. Reas¬ 
suringly, however, the other simulations soon “catch up” as their 
outer-disk star formation proceeds and more material is entrained, 
until the long-timescale outflow rates agree well. 
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Figure Dl. Shock tube tests designed to verify that numerical in-shock cooling does not significantly affect our results at our typical simulation resolution. 
See Appendixj^for details. From left to right, three idealized shock tube problems are shown which represent a BAL wind encountering a cold, dense ISM at 
radii /? = {1,10, lOOjpc, respectively, and evolved to r = 50kyrs. Different colors correspond to the particle masses (as labeled); the highest-resolution case is 
comparable to our simulations in the main text. At lower resolutions there is noticeable numerical shock-broadening, however the post-shock temperatures are 
still well-converged (i.e. they are not affected by in-shock cooling, which would systematically change the post-shock temperatures at different resolutions). 


APPENDIX D: RESOLVING IN-SHOCK COOLING: 
NUMERICAL & RESOLUTION REQUIREMENTS 

The SPH hydro solver employed in this work captures the jump 
conditions associated with shocks over a finite width of order the 
kernel softening length. Gas particles passing through a numeri¬ 
cally broadened shock can radiate away energy through traditional 
cooling channels. The post-shock gas temperature can therefore be 
numerically reduced via this “in-shock cooling" effect when shocks 
are broadened (e.g., [Hutchings & Thomas ||2000] jCreasey et al.j 
[201 Ij ). In-shock cooling can become significant when the cool¬ 
ing timescale for gas moving through the shock is comparable to 
the resolution-dependent shock crossing timescale. For the case of 
BAL winds with input velocities of ~ 10"^ km s“\ the post-shock 
gas is expected to be heated to of order T ~ 10^ K where it will cool 
inefficiently owning to two-temperature plasma effects ( jFaucher-j 
[Giguere & Quataert|2012j l. In the absence of efficient cooling chan¬ 
nels for the post-shock gas, the outfiows will remain energy con¬ 
serving, efficiently driving outfiows via PdV work on the ambient 
ISM. As a result, any significant amount of in-shock cooling can 
impact the post-shock gas temperature, and thus numerically mod¬ 
ulate the quasar feedback efficiency studied in this paper. In this 
section, we investigate the magnitude of in-shock cooling effects 
via idealized numerical experiments and find that in-shock cool¬ 
ing should be minimal for the appropriate physical conditions and 
resolutions used throughout this paper. 

In the full feedback simulations presented in this paper, BAL 
winds are implemented by imparting kicks of 5,000 km s“^ to par¬ 
ticles near the central black hole based on the accretion rate. The 
wind particles shock when they reach the ambient static ISM, ther- 
malizing their kinetic energy, and giving rise to a physical situation 


^ To ensure the simulations are otherwise exactly identical, we have had to 
re-run the v5000 simulation in the “pressure-entropy” case with a number of 
small modifications to the algorithm, and on an identical node configuration 
with pre-set values for certain random number calls. This is done for all tests 
in this section. For convenience we run the test cases at 1/8 the particle 
number in the text. 


similar to that shown in Figure 1 of |Faucher-Giguere & Quataert] 
( |20I2j ). To properly capture the full impact of BAL wind injec¬ 
tion on quasar outfiows, it is important that the thermalization of 
the BAL wind kinetic energy at the reverse shock does not suf¬ 
fer substantial from in-shock cooling. Gas densities, temperatures, 
and velocities for the reverse shock are set by the pre-shock BAL 
wind material, which is assumed to be free streaming. The non- 
homogenous density structure of the ISM and variability of the 
AGN radiation field make identifying the impact of in-shock cool¬ 
ing difficult in the full simulations directly. We instead construct 
idealized shock tube tests to recreate these conditions in a setting 
where resolution-dependent in-shock cooling can be directly iden¬ 
tified. 

We use a three-dimensional shock tube to explore the reverse 
shock density and temperature profile as a function of physical con¬ 
ditions (i.e. pre-shock density) and numerical resolution. The ide¬ 
alized initial conditions for the reverse shock include a fast moving 
medium (imitating the BAL wind material) moving into a static 
medium (imitating the ambient ISM). The BAL wind material is 
given an initial temperature of Tbal = 10"^ K, however this can 
change rapidly at the onset of the simulation due to Compton heat¬ 
ing off of the AGN radiation if the gas density is sufficiently low. 
The BAL material is given a velocity of vbal = 5,000 km s“^ The 
ambient ISM material is given an initial temperature of Tism = 10^ 
K and is initially static. The initial density for the BAL wind, initial 
density for the ambient ISM, and incident flux of AGN radiation 
are dependent on the location of the shock. We approximate the 
density of the pre-shock free-streaming BAL wind material to be 
given by 

(sOOOk^^) 

the density of the ambient ISM as 
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and the incident AGN radiation flux as 



erg ( L 



R 

1 pc 


-2 


s cm^ \ 10"^^ erg s~^ 


We assume flducial values ofM — IMo/yr, v = 5000 km/sec, and 
L — 10"^^ erg/sec and run tests for R — {1,10,100} pc. The shock 
tube uses periodic boundary conditions in a rectangular prism of 
dimension I xR/ (800 pc) xR/ (800 pc) kpc. 

Figure |dT] shows the gas density and temperature profiles 
across the idealized shock at t = 50 kyrs. The three panels show 
different values for the ambient gas density and AGN radiation 
flux (corresponding to /? = {1,10,100} pc, as described in the pre¬ 
vious paragraph) tests with the legend indicating the gas particle 
mass resolution in each test. We And that the lowest resolution test 
(black line) blurs the location of the reverse shock substantially in 
the /? = 10 and R = 100 pc tests. The two higher resolution tests 
present with less blurring the of the reverse shock, however there 
is still an offset present in the location of the reverse shock ow¬ 
ing to the low particle number in the pre shock low density BAL 
wind material. In terms of in-shock cooling, we And that the post 
shock gas forms a stable and nearly flat temperature profile which 
shows little variation as we change the mass/particle resolution for 
our highest two resolution tests. Although some shock broadening 
is present, the post-shock gas temperatures are not strongly (if at 
all) impacted by the increasing resolution. If in-shock cooling were 
present, we would instead expect the post-shock gas temperature 
to decrease with lower mass resolution. Since the high resolution 
tests explored here (m = lOM©) have resolution comparable to that 
used in the simulations presented in this paper and show little indi¬ 
cation in-shock cooling, we conclude that in-shock cooling should 
not significantly impact the post-shock gas temperatures, and there¬ 
for should not have a significant impact on the results presented in 
this paper. 
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